I’m looking forward to taking part in Introduction to Open Data Science course. I hope to learn how to use open and free software like R to analyze data as well as to learn more about different statistical methods.
# This is a so-called "R chunk" where you can write R code.
date()
## [1] "Thu Nov 19 18:34:22 2020"
date()
## [1] "Thu Nov 19 18:34:22 2020"
lrn14<-read.table("http://s3.amazonaws.com/assets.datacamp.com/production/course_2218/datasets/learning2014.txt", , sep=",", header=TRUE)
dim(lrn14)
## [1] 166 7
str(lrn14)
## 'data.frame': 166 obs. of 7 variables:
## $ gender : chr "F" "M" "F" "M" ...
## $ age : int 53 55 49 53 49 38 50 37 37 42 ...
## $ attitude: num 3.7 3.1 2.5 3.5 3.7 3.8 3.5 2.9 3.8 2.1 ...
## $ deep : num 3.58 2.92 3.5 3.5 3.67 ...
## $ stra : num 3.38 2.75 3.62 3.12 3.62 ...
## $ surf : num 2.58 3.17 2.25 2.25 2.83 ...
## $ points : int 25 12 24 10 22 21 21 31 24 26 ...
The dataset lerning2014 includes 7 different variables and 166 observations. Observations represent students. Besides background variables gender and age, the dataset includes information how student scored in a exam (variable points) and their attidude towards statistics on a scale 1 to 5 (attitude). The variable deep is the average points on a scale 1 to 5 from questions concerning deep learning approaches. Variables stra and surf include information on average points from questions about strategic learning and surface learning approaches.
More information about the original dataset, that was used to form this dataset, can be found here
library(GGally)
## Loading required package: ggplot2
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(ggplot2)
p <- ggpairs(lrn14, mapping = aes(col =gender, alpha = 0.3), lower = list(combo = wrap("facethist", bins = 20)))
p
summary(lrn14)
## gender age attitude deep
## Length:166 Min. :17.00 Min. :1.400 Min. :1.583
## Class :character 1st Qu.:21.00 1st Qu.:2.600 1st Qu.:3.333
## Mode :character Median :22.00 Median :3.200 Median :3.667
## Mean :25.51 Mean :3.143 Mean :3.680
## 3rd Qu.:27.00 3rd Qu.:3.700 3rd Qu.:4.083
## Max. :55.00 Max. :5.000 Max. :4.917
## stra surf points
## Min. :1.250 Min. :1.583 Min. : 7.00
## 1st Qu.:2.625 1st Qu.:2.417 1st Qu.:19.00
## Median :3.188 Median :2.833 Median :23.00
## Mean :3.121 Mean :2.787 Mean :22.72
## 3rd Qu.:3.625 3rd Qu.:3.167 3rd Qu.:27.75
## Max. :5.000 Max. :4.333 Max. :33.00
From the graph above we can see the distrubutions of different variables by gender. Below the graph there are summary statistics for different variables. The avereage age of study participants is 25.51 years. Majority of participants are women. Attitude variable has the strongest correlation with exam points and the scatter plot of these two variables also indicate that there may be some kind assocation between these variables.
model1 <- lm(points ~ attitude + stra + surf, data = lrn14)
summary(model1)
##
## Call:
## lm(formula = points ~ attitude + stra + surf, data = lrn14)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.1550 -3.4346 0.5156 3.6401 10.8952
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.0171 3.6837 2.991 0.00322 **
## attitude 3.3952 0.5741 5.913 1.93e-08 ***
## stra 0.8531 0.5416 1.575 0.11716
## surf -0.5861 0.8014 -0.731 0.46563
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.296 on 162 degrees of freedom
## Multiple R-squared: 0.2074, Adjusted R-squared: 0.1927
## F-statistic: 14.13 on 3 and 162 DF, p-value: 3.156e-08
model2 <- lm(points ~ attitude , data = lrn14)
summary(model2)
##
## Call:
## lm(formula = points ~ attitude, data = lrn14)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.9763 -3.2119 0.4339 4.1534 10.6645
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.6372 1.8303 6.358 1.95e-09 ***
## attitude 3.5255 0.5674 6.214 4.12e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.32 on 164 degrees of freedom
## Multiple R-squared: 0.1906, Adjusted R-squared: 0.1856
## F-statistic: 38.61 on 1 and 164 DF, p-value: 4.119e-09
We choose variables _attitude, stra and surf as independents variables to fit a linear regression model where points i.e. exam points is the dependent variable. The chosen independent variables were the most correlated variables with the dependent variable. From Model 1 summary we can see that only attitude has a statistically significant coefficient (3.3952) as its p-value <0.0001. With low p-value we can reject the null hypothesis that the coefficent is zero. When we fit Model 2 without stra and surf the coefficent for attitude is still significant and its value (3.5255) is close to that in Model 1.
The multiple R-squared tells us how much our models explain the variation in points. Multiple R-squared are quite close to each other in both models (Model 1 0.2074 and Model 2 0.1906 ) so roughly both models explain approx. 20% of the varioation in points. In Model 2 with a single independent variable the miltiple R-squared is basically the square of correlation coefficient between the independent and the dependent variable.
par(mfrow =c(2,2))
plot(model2, which = c(1,2,5) )
From the QQ plot we can see that the residuals reasonable well fall in with the line which would indicate that the assumption of normality of error terms are not violated.
From the plot where residuals and fitted values from the model are plotted we can see that there is no clear pattern. The size of error term does not depend on the independent variables which indicates that the constant variance of error term assumption is not violated.
The leverage plot helps us to estimate if there are influential outliers in our data that affect out estimation results. From our leverage plot we see that there are no influential outliers, the x-axis scale is very low.
date()
## [1] "Thu Nov 19 18:34:29 2020"
library(dplyr); library(ggplot2)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
setwd("~/R/IODS-project")
pormat<-read.csv("data/pormat.csv", header=T)
str(pormat)
## 'data.frame': 370 obs. of 35 variables:
## $ school : chr "GP" "GP" "GP" "GP" ...
## $ sex : chr "F" "F" "F" "F" ...
## $ age : int 15 15 15 15 15 15 15 15 15 15 ...
## $ address : chr "R" "R" "R" "R" ...
## $ famsize : chr "GT3" "GT3" "GT3" "GT3" ...
## $ Pstatus : chr "T" "T" "T" "T" ...
## $ Medu : int 1 1 2 2 3 3 3 2 3 3 ...
## $ Fedu : int 1 1 2 4 3 4 4 2 1 3 ...
## $ Mjob : chr "at_home" "other" "at_home" "services" ...
## $ Fjob : chr "other" "other" "other" "health" ...
## $ reason : chr "home" "reputation" "reputation" "course" ...
## $ guardian : chr "mother" "mother" "mother" "mother" ...
## $ traveltime: int 2 1 1 1 2 1 2 2 2 1 ...
## $ studytime : int 4 2 1 3 3 3 3 2 4 4 ...
## $ schoolsup : chr "yes" "yes" "yes" "yes" ...
## $ famsup : chr "yes" "yes" "yes" "yes" ...
## $ activities: chr "yes" "no" "yes" "yes" ...
## $ nursery : chr "yes" "no" "yes" "yes" ...
## $ higher : chr "yes" "yes" "yes" "yes" ...
## $ internet : chr "yes" "yes" "no" "yes" ...
## $ romantic : chr "no" "yes" "no" "no" ...
## $ famrel : int 3 3 4 4 4 4 4 4 4 4 ...
## $ freetime : int 1 3 3 3 2 3 2 1 4 3 ...
## $ goout : int 2 4 1 2 1 2 2 3 2 3 ...
## $ Dalc : int 1 2 1 1 2 1 2 1 2 1 ...
## $ Walc : int 1 4 1 1 3 1 2 3 3 1 ...
## $ health : int 1 5 2 5 3 5 5 4 3 4 ...
## $ failures : int 0 1 0 0 1 0 1 0 0 0 ...
## $ paid : chr "yes" "no" "no" "no" ...
## $ absences : int 3 2 8 2 5 2 0 1 9 10 ...
## $ G1 : int 10 10 14 10 12 12 11 10 16 10 ...
## $ G2 : int 12 8 13 10 12 12 6 10 16 10 ...
## $ G3 : int 12 8 12 9 12 12 6 10 16 10 ...
## $ alc_use : num 1 3 1 1 2.5 1 2 2 2.5 1 ...
## $ high_use : logi FALSE TRUE FALSE FALSE TRUE FALSE ...
The data set pormat includes 35 different variables and 370 observations. The data set has been constructed by joining two data sets that include information on students from two Portuguese schools. The two original data sets included information on students performance in mathematics and Portuguese language. Variables G1-G3 are the grades and when data sets were joined the averages of grades of these two subjects were calculated to create grades variables in the joined data set. Also averages of the numbers absences and failures from these two subject were calculated when forming the joined data set. paid variable contains information if the student got paid extra classes in the subject, in the joined data set paid variable represent the answer to the question if the student got paid classes in Portuguese language. Other variables are background variables such as age, sex, parents’ level of education etc.
When the joined data set was constructed, two new variables were created alc_use which is the average of workday (Dalc) and weekend (Walc) alcohol consumption and high_use which is binary variable describing if the average alcohol consumption was higher than 2.
The two original datasets and the detail description of all variables in the original data sets are available here.
Absences: Students who have a lot of absences from class may also have a high alcohol consumption as absences may reflect the fact that student skip school because they feel too hangover.
Failures: A high alcohol consumption may harm school performance and student who consume high amount alcohol may not pass exams as often as other students.
Sex: It is a quite common known fact that men consume more alcohol than women.
PStatus (parent’s cohabitation status): Students whose parents live apart may consume more alcohol because their single parents possibly have to work a lot to provide for the family and can’t spent time to supervise what children are doing or a high alcohol consumption may also be a child’s reaction to parents’ separation.
library(tidyr); library(dplyr); library(ggplot2); library(gmodels)
vars<-subset(pormat, select= c("failures","sex","absences", "Pstatus", "alc_use", "high_use"))
# draw a bar plot of each variable
gather(vars) %>% ggplot(aes(value)) + facet_wrap("key", scales = "free")+geom_bar()
g1<- ggplot(pormat, aes(x=high_use, y=absences,))
g1 + geom_boxplot() + ylab("absences")+ggtitle("Student absences by alcohol consumption and sex")
table(high_use = pormat$high_use, failures = pormat$failures)
## failures
## high_use 0 1 2 3
## FALSE 238 12 8 1
## TRUE 87 12 9 3
CrossTable(pormat$high_use, pormat$sex, prop.c=TRUE, prop.r=F, prop.chisq=F, prop.t=F)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Col Total |
## |-------------------------|
##
##
## Total Observations in Table: 370
##
##
## | pormat$sex
## pormat$high_use | F | M | Row Total |
## ----------------|-----------|-----------|-----------|
## FALSE | 154 | 105 | 259 |
## | 0.790 | 0.600 | |
## ----------------|-----------|-----------|-----------|
## TRUE | 41 | 70 | 111 |
## | 0.210 | 0.400 | |
## ----------------|-----------|-----------|-----------|
## Column Total | 195 | 175 | 370 |
## | 0.527 | 0.473 | |
## ----------------|-----------|-----------|-----------|
##
##
CrossTable(pormat$high_use, pormat$Pstatus, prop.c=TRUE, prop.r=F, prop.chisq=F, prop.t=F)
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Col Total |
## |-------------------------|
##
##
## Total Observations in Table: 370
##
##
## | pormat$Pstatus
## pormat$high_use | A | T | Row Total |
## ----------------|-----------|-----------|-----------|
## FALSE | 26 | 233 | 259 |
## | 0.684 | 0.702 | |
## ----------------|-----------|-----------|-----------|
## TRUE | 12 | 99 | 111 |
## | 0.316 | 0.298 | |
## ----------------|-----------|-----------|-----------|
## Column Total | 38 | 332 | 370 |
## | 0.103 | 0.897 | |
## ----------------|-----------|-----------|-----------|
##
##
From the first plot we can see the distributions of our chosen variables. For example, most students’ parents live together and there are more women in our sample than men.
From the second graph we see that the median number of absences is higher among those who consume high amounts of alcohol and also limits for IQR are higher.
From the first table we see that there are more high consumers among those who have several failures in exams.
From the second table we see that 40% of men and 21% of women are have high alcohol consumption.
From the third table we see that 31% of students who parents live apart and 30% of students who parents live together are high consumers.
m1 <- glm(high_use ~ failures + absences + sex + Pstatus, data = pormat, family = "binomial")
# print out a summary of the model
summary(m1)
##
## Call:
## glm(formula = high_use ~ failures + absences + sex + Pstatus,
## family = "binomial", data = pormat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1317 -0.8455 -0.5912 1.0272 2.0345
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.02222 0.43765 -4.621 3.83e-06 ***
## failures 0.59686 0.20720 2.881 0.00397 **
## absences 0.09301 0.02336 3.982 6.84e-05 ***
## sexM 0.99666 0.24726 4.031 5.56e-05 ***
## PstatusT 0.08764 0.40210 0.218 0.82746
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 452.04 on 369 degrees of freedom
## Residual deviance: 406.95 on 365 degrees of freedom
## AIC: 416.95
##
## Number of Fisher Scoring iterations: 4
# print out the coefficients of the model
coef(m1)
## (Intercept) failures absences sexM PstatusT
## -2.02221996 0.59686177 0.09301028 0.99666359 0.08764120
# compute odds ratios (OR)
OR <- coef(m1) %>% exp
# compute confidence intervals (CI)
CI<-confint(m1)%>%exp
## Waiting for profiling to be done...
# print out the odds ratios with their confidence intervals
cbind(OR, CI)
## OR 2.5 % 97.5 %
## (Intercept) 0.1323613 0.05365585 0.3015602
## failures 1.8164095 1.21846564 2.7632268
## absences 1.0974730 1.05068550 1.1517051
## sexM 2.7092276 1.67857047 4.4338543
## PstatusT 1.0915964 0.50797217 2.4883339
From the summary of our model we can see that failures, absences, and sex variables have a statistically significant coefficients (p<0.001) but Psatus variable do not have a significant coefficient. From odds ratios we can see that men have 2.71 or according to confidence intervals from 1.68 to 4.43 times higher odds of high alcohol consumption than women. Number of failures and absences increase the odds of high alcohol consumption by 81.6% and 9.7% respectively for every one unit increase in the value of the variable. From the confidence interval of Pstatus we can make the same conclusion as from the p-value. The CI includes the value one which indicates that the OR is not statistically significant. To conclude failures, absences and sex seem to be associated with a high alcohol consumption as hypothezized previously but parent’s cohabitation status is not associated with a high alcohol consumption.
m2 <- glm(high_use ~ failures + absences + sex , data = pormat, family = "binomial")
# predict() the probability of high_use
probabilities <- predict(m2, type = "response")
# add the predicted probabilities to data
pormat <- mutate(pormat, probability = probabilities)
# use the probabilities to make a prediction of high_use
pormat <- mutate(pormat, prediction = probability>0.5)
# tabulate the target variable versus the predictions
table(high_use = pormat$high_use, prediction = pormat$prediction)
## prediction
## high_use FALSE TRUE
## FALSE 252 7
## TRUE 78 33
table(high_use = pormat$high_use, prediction = pormat$prediction)%>%prop.table()%>%addmargins
## prediction
## high_use FALSE TRUE Sum
## FALSE 0.68108108 0.01891892 0.70000000
## TRUE 0.21081081 0.08918919 0.30000000
## Sum 0.89189189 0.10810811 1.00000000
# plot the target variable versus the predictions
g <- ggplot(pormat, aes(x = high_use, y = probability))
g+geom_point()
loss_func <- function(class, prob) {
n_wrong <- abs(class - prob) > 0.5
mean(n_wrong)
}
# call loss_func to compute the average number of wrong predictions in the (training) data
loss_func(class = pormat$high_use, prob = pormat$probability)
## [1] 0.2297297
As Pstatus was not statistically significant we leave it out from the final model. Our final model predicted 252 FALSE cases out of 259 right as for TRUE cases only 33 out of 111 was predicted right. The total proportion of inaccurately classified individuals was 23%.
library(boot)
cv <- cv.glm(data = pormat, cost = loss_func, glmfit = m2, K = 10)
# average number of wrong predictions in the cross validation
cv$delta[1]
## [1] 0.2486486
According the 10-fold cross validation out model has about 0.24 error. That is a bit smaller than in the DataCamp model although our final model is exactly the same. However our data has fewer observations so the comparisons is not so straightforward.
date()
## [1] "Thu Nov 19 18:34:31 2020"
library(tidyr); library(dplyr); library(ggplot2); library(gmodels)
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
# load the data
data("Boston")
# explore the dataset
dim(Boston)
## [1] 506 14
str(Boston)
## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
The Boston data set has 506 observations and 14 variables. The data include information on housing values in suburbs of Boston. Full list of variables and descriptions of variables are available here
pairs(Boston)
pairs(Boston[c("crim","age","dis","rad","tax","black","medv")])
summary(Boston)
## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00
From the graphs we can see the distributions of variables and their relationship with other variables. When we take a closer look some selected variables, we see that for example crime rates are higher in areas where there are larger proportions of older houses or in areas where the mean of distances to employment centers is shorter.
# center and standardize variables
boston_scaled <- scale(Boston)
# summaries of the scaled variables
summary(boston_scaled)
## crim zn indus chas
## Min. :-0.419367 Min. :-0.48724 Min. :-1.5563 Min. :-0.2723
## 1st Qu.:-0.410563 1st Qu.:-0.48724 1st Qu.:-0.8668 1st Qu.:-0.2723
## Median :-0.390280 Median :-0.48724 Median :-0.2109 Median :-0.2723
## Mean : 0.000000 Mean : 0.00000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.007389 3rd Qu.: 0.04872 3rd Qu.: 1.0150 3rd Qu.:-0.2723
## Max. : 9.924110 Max. : 3.80047 Max. : 2.4202 Max. : 3.6648
## nox rm age dis
## Min. :-1.4644 Min. :-3.8764 Min. :-2.3331 Min. :-1.2658
## 1st Qu.:-0.9121 1st Qu.:-0.5681 1st Qu.:-0.8366 1st Qu.:-0.8049
## Median :-0.1441 Median :-0.1084 Median : 0.3171 Median :-0.2790
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.5981 3rd Qu.: 0.4823 3rd Qu.: 0.9059 3rd Qu.: 0.6617
## Max. : 2.7296 Max. : 3.5515 Max. : 1.1164 Max. : 3.9566
## rad tax ptratio black
## Min. :-0.9819 Min. :-1.3127 Min. :-2.7047 Min. :-3.9033
## 1st Qu.:-0.6373 1st Qu.:-0.7668 1st Qu.:-0.4876 1st Qu.: 0.2049
## Median :-0.5225 Median :-0.4642 Median : 0.2746 Median : 0.3808
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 1.6596 3rd Qu.: 1.5294 3rd Qu.: 0.8058 3rd Qu.: 0.4332
## Max. : 1.6596 Max. : 1.7964 Max. : 1.6372 Max. : 0.4406
## lstat medv
## Min. :-1.5296 Min. :-1.9063
## 1st Qu.:-0.7986 1st Qu.:-0.5989
## Median :-0.1811 Median :-0.1449
## Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.6024 3rd Qu.: 0.2683
## Max. : 3.5453 Max. : 2.9865
# change the object to data frame
boston_scaled<-as.data.frame(boston_scaled)
# create a quantile vector of crim
bins <- quantile(boston_scaled$crim)
# create a categorical variable 'crime'
crime <- cut(boston_scaled$crim, breaks = bins, include.lowest = TRUE, label=c("low","med_low","med_high","high"))
# look at the table of the new factor crime
table(crime)
## crime
## low med_low med_high high
## 127 126 126 127
# add the new categorical value to scaled data
boston_scaled <- data.frame(boston_scaled, crime)
# remove the crime variable from the data set
boston_scaled <- dplyr::select(boston_scaled, -crim)
# number of rows in the Boston dataset
n <- nrow(boston_scaled)
# choose randomly 80% of the rows
ind <- sample(n, size = n * 0.8)
# create train set
train <- boston_scaled[ind,]
# create test set
test <- boston_scaled[-ind,]
After standardization, we can see from the summary table that all variables now have the mean zero and the minimum and maximum values of all scaled variables varies in much smaller intervals than in the original data.
After creating a categorical variable for crime rates, we can see from the table that observations are quite equally distributed in four categories.
# linear discriminant analysis
lda.fit <- lda(crime~., data = train)
# print the lda.fit object
lda.fit
## Call:
## lda(crime ~ ., data = train)
##
## Prior probabilities of groups:
## low med_low med_high high
## 0.2425743 0.2500000 0.2425743 0.2648515
##
## Group means:
## zn indus chas nox rm age
## low 0.9603007 -0.8877158 -0.11163110 -0.8678412 0.4322590 -0.8570834
## med_low -0.1098365 -0.2463064 -0.07742312 -0.5301939 -0.1853758 -0.2457814
## med_high -0.3804849 0.1637265 0.08924136 0.4063825 -0.0191873 0.4237528
## high -0.4872402 1.0170108 -0.05155709 1.0770807 -0.3833050 0.8059241
## dis rad tax ptratio black lstat
## low 0.8302295 -0.6865511 -0.7608232 -0.444665440 0.38499661 -0.77516658
## med_low 0.3601183 -0.5486376 -0.4587498 0.005447313 0.32135485 -0.07266494
## med_high -0.3449008 -0.4111534 -0.3038295 -0.288654539 0.07858534 0.08529259
## high -0.8500432 1.6392096 1.5148289 0.782035629 -0.82427597 0.88654355
## medv
## low 0.52854368
## med_low -0.06266871
## med_high 0.08685629
## high -0.71508722
##
## Coefficients of linear discriminants:
## LD1 LD2 LD3
## zn 0.094109369 0.90085820 -0.81734309
## indus -0.035554476 -0.15920051 0.32687687
## chas 0.025747812 -0.09912745 0.03586020
## nox 0.459379374 -0.63925638 -1.37302914
## rm 0.049755854 0.01301499 -0.13551227
## age 0.247839076 -0.34734052 -0.04417047
## dis -0.081490631 -0.46593622 0.32020903
## rad 3.157275451 0.98636087 -0.04463778
## tax 0.005525413 -0.21033169 0.59500414
## ptratio 0.153101786 0.07396940 -0.16671214
## black -0.129431041 0.02863432 0.14176732
## lstat 0.194440033 -0.29804647 0.37712029
## medv 0.065744528 -0.41702050 -0.15389035
##
## Proportion of trace:
## LD1 LD2 LD3
## 0.9526 0.0354 0.0120
# target classes as numeric
classes <- as.numeric(train$crime)
#plot LDA biplot
plot(lda.fit, dimen = 2, col = classes, pch = classes)
From the output of our LDA model, we see that the first linear discriminant explain 94.5% of between group variance. This can be also seen from the plot where along the x axis (LD1) the separation of different groups is somewhat clearer, especially for the high category, than along the y axis (LD2).
# save the correct classes from test data
correct_classes <- test$crime
# remove the crime variable from test data
test <- dplyr::select(test, -crime)
# predict classes with test data
lda.pred <- predict(lda.fit, newdata = test)
# cross tabulate the results
table(correct = correct_classes, predicted = lda.pred$class)
## predicted
## correct low med_low med_high high
## low 17 12 0 0
## med_low 7 16 2 0
## med_high 1 6 20 1
## high 0 0 0 20
From the table we see that our LDA model did fairly good job as the majority of predictions were correct in most categories. The med low category was the most hardest to predict right.
set.seed(123)
library(MASS); library(ggplot2)
data("Boston")
# center and standardize variables
boston_scaled <- scale(Boston)
boston_scaled<-as.data.frame(boston_scaled)
# euclidean distance matrix
dist_eu <- dist(boston_scaled)
# look at the summary of the distances
summary(dist_eu)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1343 3.4625 4.8241 4.9111 6.1863 14.3970
# k-means clustering:
km <-kmeans(boston_scaled, centers = 2)
pairs(boston_scaled, col = km$cluster)
pairs(boston_scaled[6:10], col = km$cluster)
# Investigate the optimal number of clusters with the total of within cluster sum of squares (tWCSS)
# determine the number of clusters
k_max <- 10
# calculate the total within sum of squares
twcss <- sapply(1:k_max, function(k){kmeans(boston_scaled, k)$tot.withinss})
# visualize the results
qplot(x = 1:k_max, y = twcss, geom = 'line')
# k-means clustering: 3 clusters
km <-kmeans(boston_scaled, centers = 3)
pairs(boston_scaled, col = km$cluster)
pairs(boston_scaled[6:10], col = km$cluster)
When we first run K-means clustering with two clusters we can see for example rad and tax variables seem have a clear effect on the clustering results. When we investigate the optimal nubmer of clusters, three clusters seem to be enough according the total of within cluster sum of squares. With three clusters again rad and tax seem to affect the clustering results.
library(MASS)
data("Boston")
set.seed(123)
# center and standardize variables
boston_scaled <- scale(Boston)
boston_scaled<-as.data.frame(boston_scaled)
km <-kmeans(boston_scaled, centers = 3)
boston_scaled1<-data.frame(boston_scaled, km$cluster)
# linear discriminant analysis
lda.fit2 <- lda(km.cluster~., data = boston_scaled1)
# print the lda.fit object
lda.fit2
## Call:
## lda(km.cluster ~ ., data = boston_scaled1)
##
## Prior probabilities of groups:
## 1 2 3
## 0.2806324 0.3992095 0.3201581
##
## Group means:
## crim zn indus chas nox rm
## 1 0.9693718 -0.4872402 1.074440092 -0.02279455 1.04197430 -0.4146077
## 2 -0.3549295 -0.4039269 0.009294842 0.11748284 0.01531993 -0.2547135
## 3 -0.4071299 0.9307491 -0.953383032 -0.12651054 -0.93243813 0.6810272
## age dis rad tax ptratio black
## 1 0.7666895 -0.8346743 1.5010821 1.4852884 0.73584205 -0.7605477
## 2 0.3096462 -0.2267757 -0.5759279 -0.4964651 -0.09219308 0.2473725
## 3 -1.0581385 1.0143978 -0.5976310 -0.6828704 -0.53004055 0.3582008
## lstat medv
## 1 0.85963373 -0.6874933
## 2 0.09168925 -0.1052456
## 3 -0.86783467 0.7338497
##
## Coefficients of linear discriminants:
## LD1 LD2
## crim 0.03654114 0.20373943
## zn -0.08346821 0.34784463
## indus -0.32262409 -0.12105014
## chas -0.04761479 -0.13327215
## nox -0.13026254 0.15610984
## rm 0.13267423 0.44058946
## age -0.11936644 -0.84880847
## dis 0.23454618 0.58819732
## rad -1.96894437 0.57933028
## tax -1.10861600 0.53984421
## ptratio -0.13087741 -0.02004405
## black 0.15432491 -0.06106305
## lstat -0.14002173 0.14786473
## medv 0.02559139 0.37307811
##
## Proportion of trace:
## LD1 LD2
## 0.8999 0.1001
# target classes as numeric
classes <- as.numeric(boston_scaled1$km.cluster)
# the function for lda biplot arrows
lda.arrows <- function(x, myscale = 1, arrow_heads = 0.1, color = "red", tex = 0.75, choices = c(1,2)){
heads <- coef(x)
arrows(x0 = 0, y0 = 0,
x1 = myscale * heads[,choices[1]],
y1 = myscale * heads[,choices[2]], col=color, length = arrow_heads)
text(myscale * heads[,choices], labels = row.names(heads),
cex = tex, col=color, pos=3)
}
#plot LDA biplot
plot(lda.fit2, dimen = 2, col = classes, pch = classes)
lda.arrows(lda.fit2, myscale = 3)
From the plot we can see that rad and tax varaiables are mostly correlated with the first linear discriminant and dis and age with the second linear discriminant.
set.seed(123)
model_predictors <- dplyr::select(train, -crime)
# check the dimensions
dim(model_predictors)
## [1] 404 13
dim(lda.fit$scaling)
## [1] 13 3
# matrix multiplication
matrix_product <- as.matrix(model_predictors) %*% lda.fit$scaling
matrix_product <- as.data.frame(matrix_product)
#Next, install and access the plotly package. Create a 3D plot (Cool!) of the columns of the matrix product by typing the code below.
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
classes <- as.numeric(train$crime)
plot_ly(x = matrix_product$LD1, y = matrix_product$LD2, z = matrix_product$LD3, type= 'scatter3d', mode='markers', color=classes)
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
I didn’t figure out anymore how to add clusters as colors.